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We study the evolution of linear cosmological perturbations in f(R) models of accelerated ex- 
pansion in the physical frame where the gravitational dynamics are fourth order and the matter is 
minimally coupled. These models predict a rich and testable set of linear phenomena. For each 
expansion history, fixed empirically by cosmological distance measures, there exists two branches 
of f(R) solutions that are parameterized by B oc d 2 f/dR 2 . For B < 0, which include most of the 
models previously considered, there is a short-timescale instability at high curvature that spoils 
agreement with high redshift cosmological observables. For the stable B > branch, f(R) models 
can reduce the large-angle CMB anisotropy, alter the shape of the linear matter power spectrum, 
and qualitatively change the correlations between the CMB and galaxy surveys. All of these phe- 
nomena are accessible with current and future data and provide stringent tests of general relativity 
on cosmological scales. 

PACS numbers: 



I. INTRODUCTION 

Cosmic acceleration can be explained either by miss- 
ing energy with an exotic equation of state, dubbed dark 
energy, or by a modification of gravity on large scales. 
Indeed the cosmological constant can be considered ei- 
ther as a constant added to the Einstein-Hilbert action 
or as vacuum energy. Non-trivial modifications where the 
addition is a non- linear function /(i?) of the Ricci scalar 
that becomes important only at the cosmologically low 
values of R have also been shown to cause acceleration 
0, 0, 3 • They are furthermore free of ghosts and other 
types of instabilities for a wide range of interesting cases 

am. 

Solar-system tests of gravity provide what is perhaps 
the leading challenge to f(R) models as a complete the- 
ory of gravity @. The equivalence of f(R) models to 
scalar-tensor theories lead to conflicts with parameter- 
ized post-Newtonian constraints at a background cosmo- 
logical density of matter. It is however still controversial 
whether the whole class of f(R) modifications can be 
ruled out by this equivalence. Matter in the solar sys- 
tem becomes non-minimally coupled in the transformed 
frame leading to non-trivial modifications of the scalar 
field potential. In the original Jordan — or physical — 
frame, it has been shown that the Schwarzschild metric 
solves the modified Einstein equations of a wide range 
of f(R) models but this solution is not necessarily 

relevant for the solar system 0. Recent work has also 
raised the question as to whether solar-system gravity 
problems may become tractable if f(R) is viewed as sim- 
ply a first-order correction term to the high R limit of 
general relativity [H [Q, 0, 03 • 



Regardless of the outcome of small-scale tests of grav- 
ity in f(R) models, it is worthwhile to examine the cos- 
mological consequences of treating f(R) as an effective 
theory valid for a cosmologically appropriate range of 
curvatures. At the very least by making concrete pre- 
dictions of cosmological phenomena in these models, one 
gains insight on how cosmology can test gravity at the 
largest scales. 

In this Paper, we develop linear perturbation theory 
for predicting cosmological observables such as the Cos- 
mic Microwave Background (CMB) and the large-scale 
structure of the universe exhibited in galaxy surveys. We 
work in the physical frame where the matter is minimally 
coupled and obeys simple conservation laws. 

We begin in S|TI] by reviewing the properties of f(R) 
models and their relationship to the expansion history of 
the universe. In ijllll we derive the fourth order pertur- 
bation equations and, using general properties demanded 
by energy- momentum conservation Il5l . Il6j . recast 
them into a tractable second order form. In SjlVI we iden- 
tify a short time scale instability that renders a wide class 
of f(R) models not viable cosmologically. We present so- 
lutions on the stable branch in SjV] and explore their im- 
pact on cosmological power spectra in WI1 We discuss 
these results in Will 



II. EXPANSION HISTORY 

We consider a modification to the Einstein-Hilbert ac- 
tion of the form [ijj 
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where R is the Ricci scalar, which we will sometimes refer 
to as the curvature, ^ 2 = 8nG, and C m is the matter 
Lagrangian. Variation with respect to the metric yields 
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the modified Einstein equations 



G a p + fuRap - ( 2 - OfR)g a p - V„V^/fl = [i T a p,(2) 

where f R = df/dR and likewise f RR = d 2 f/dR 2 below. 
We define the metric to include scalar linear perturba- 
tions around a flat FRW background in the Newtonian 
or longitudinal gauge 



ds 2 = -(1 + 2V)dt 2 + a 2 {l + 2<i>)da 



(3) 



Given that the expansion history and dynamics of linear 
perturbations are well-tested in the high curvature, high 
redshift limit by the CMB, we restrict our considerations 
to models that satisfy [l], 0, HJ 



lim f(R)/R-+0. 



(4) 



With this restriction, the main modifications for viable 
models with stable high curvature limits arise well af- 
ter the radiation becomes a negligible contributor to the 
stress energy tensor. We can then take it to have the 
matter-dominated form 
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The modified Einstein equations with the FRW back- 
ground metric and 5 = q = yields the modified Fried- 
mann equation 



H 2 f R (HH' + H 2 ) + \f + H 2 f RR R> = 

6 3 



(G) 



Here and throughout, primes denote derivatives with re- 
spect to In a. 

There is sufficient freedom in the function f(R) to re- 
produce any desired expansion history H. Hence the ex- 
pansion history alone cannot be used as a test of general 
relativity though it can rule out specific forms of f{R). 
The dynamics of linear perturbations on the other hand 
do test general relativity as we shall see. 

We therefore seek to determine a family of f(R) func- 
tions that is consistent with a given expansion history 
PJ, EllH HI HI- Without loss of generality, we can 
parameterize the expansion history in terms of an equiv- 
alent dark energy model 



# 2 = y(p + PD E ) 



(7) 



This yields a second order differential equation for f(R) 



f R (HH' + H 2 ) + -f + H 2 f RR R' 



At 2 yODE 
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where H 2 and R are fixed functions of In a given the 
matching to the dark energy model. 



For convenience, let us define the dimensionless quan- 
tities 



E 



^=3(AE + E>), y=J^ 



H, 



(9) 



Here H = H(\na = 0) = /i/ 2 997.9 Mpc" 1 is the Hub- 
ble constant. The modified Friedmann equation can be 
recast into an inhomogeneous differential equation for 
y(lna) 



L[y] 



p DE fAE' + E" 
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where the differential operator on the lhs is given by 

4E" + E" n 
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For illustrative purposes, we take an expansion history 
that matches a dark energy model with a constant equa- 
tion of state w, 



E = (i-n BE ) a - 3 + n BE a- 3 ( 1+ ^ 



and 



M 2 Pde 
3 H 2 



-3(1+m>) 
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Since Eq. (fTl]) is a second-order differential equation, the 
expansion history does not uniquely specify f{R) but in- 
stead allows a family of solutions that are distinguished 
by initial conditions. This additional freedom reflects the 
fourth-order nature of f(R) gravity. 

To set the initial conditions, take y± to be the two 
solutions of the homogeneous equation L[y] = 0. At high 
curvature, these solutions are power laws y± cx a p± with 



P± = 



-7±y/73 



(14) 



Since p- ~ —3.9, stimulation of this decaying mode vi- 
olates the condition that f R /R — * at high R. We 
therefore set its amplitude to zero in our solutions (c.f. 
[l^, [H, Hij]). The particular solution in the high curva- 
ture limit becomes 



en 



DE 



-3(l+io) 



(15) 



6w 2 + 5w — 2 

Therefore when numerically integrating Eq. (|11|) we take 

yQiidi) = Ay+Qncii) + y paz tQ-nai) , (16) 
y' (In en) = p + Ay + (laai) - 3(1 + w)y part (lnai) , 

at some initial epoch eti ~ 10 -2 . 

Since the modifications to gravity appear at low red- 
shifts, it is more convenient to parameterize the individ- 
ual solutions in the family by the final conditions rather 
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than the growing mode amplitude of the initial condi- 
tions. Given that a constant f(R) is simply a cosmo- 
logical constant and a linear one represents a rescaling 
of G or fi, it is Jbm, the second derivative, that controls 
phenomena that are unique to the modification. In par- 
ticular, we shall see that a specific dimensionless quantity 
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[rr p/ H 
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(17) 
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is most closely linked with the phenomenology. B is a 
strongly growing function in our solutions and in the high 
curvature limit has a growth rate 



B' 



(18) 



given by ps = 3 + p+, if the growing mode dominates, 
and pb = — 3w, if the particular mode dominates. 

We will therefore characterize solutions with a given 
expansion history family by Bq = _B(lna = 0). If Bq = 
and the background expansion is given by w = —1 
then f(R) — const., B(lna) — and the model has a 
true cosmological constant. More generally B = will 
correspond in linear theory to the dynamics of a dark 
energy component for scales above the dark energy sound 
horizon. 

In Fig.[ll we plot the family of f(R) models that match 
two representative expansion histories parameterized by 
(u>, £lr,E,h). These models are chosen to be consistent 
with current WMAP CMB data [25| and span a range 
that is consistent with supernovae acceleration measures. 
Given the similarity between these models, we will take 
the w = — 1 ACDM expansion history for illustrative 
purposes below. 

The linear perturbation analysis that follows does not 
require the matching to the specific expansion histories 
parameterized by (w, f2rjE, h) here. This reverse engi- 
neering is only a device to find observationally accept- 
able f(R) models. What is required is that the back- 
ground solutions provide _ff(lna) and _B(lna). On the 
other hand, linear perturbation theory does inform the 
choice of a background solution. We shall find in SjTV] 
that the B < branch of the family is unstable to linear 
perturbations in the high curvature regime. 



III. LINEAR PERTURBATION EQUATIONS 

The modified Einstein equations © represent a fourth- 
order set of differential equations for the two metric per- 
turbations 'J and $ in the presence of matter density and 
momentum fluctuations S and q. To solve this system of 
equations, we introduce auxiliary parameters to recast it 
as a larger set of second-order differential equations. It 
is numerically and pedagogically advantageous to choose 
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FIG. 1: Every expansion history that can be parameterized by 
a dark energy model with prjE(lna) can be reproduced by a one 
parameter family of f(R) models, indexed by So tx + fn) 

at the present epoch (left end point of curves), that approaches 
the Einstein-Hilbert action in the high curvature limit, (a) ACDM 
expansion history (to = —1, Cde = 0.76, h = 0.73). (b) Dynamical 
dark energy expansion history (w = —0.9, Hfje = 0.73, h = 0.69). 



these auxiliary parameters so that their effect vanishes 
at large scales and early times. 

On superhorizon scales (k/aH <C 1), the evolution of 
metric perturbations must be consistent with the back- 
ground evolution provided that the background solu- 
tion is valid, i.e. that fluctuations about it are stable. 
Bertschinger [l|| showed that the familiar conservation 
of the curvature fluctuation on comoving hypersurfaces 
(C' = 0) for adiabatic fluctuations in a flat universe ap- 
plies to any metric based modified gravity model that 
obeys energy momentum conservation V^T^ = 0. The 
gauge transformation into the Newtonian gauge 

C = <f> + Hq (19) 

implies 

£ = $' + H'q + Hq' = , (k = 0) , (20) 
and momentum conservation requires 

Hq' = (21) 

so that 

$'-* + if'g = 0, (fc = 0). (22) 

Combining Eq. (|2~lj) and (|22[) yields a second order differ- 
ential equation for the Newtonian metric perturbations 

El 

$"-*'-_$'-(_-_)* = 0, (fc = 0).(23) 



H W 



4 



The evolution of the metric fluctuations must in this way 
be consistent with the expansion history defined by H. 
Note that this equation applies to any modified gravi- 
tational scenario that satisfies the required conditions. 
The DGP braneworld acceleration model represents 
another valid application [27| • What does require a spec- 
ification of a theory is the relation between $ and "J. Un- 
der general relativity and assuming T^ v takes the matter- 
only form of Eq. ([5]), the closure relation is <& = — W and 
Eq. (I23p in fact applies on all scales. With a dynamical 
dark energy component in T„ v , it applies above the dark 
energy sound horizon (28j . 

To capture the metric evolution of f(R) models for 
k 0, let us introduce two parameters: 9 the deviation 
from ( conservation, Eq. (j22"|) 

H' / k \ 2 

(i = $' -ty + H'q = — — B6, (24) 



and e the deviation from the superhorizon metric evolu- 
tion Eq. (|2"5j) 



H" / N' IT" " 

<l>" _ qj> _ - ( — — I ^ 

H> Iff H' 



aH 



Be. 



(25) 



The coefficients in front of the deviation parameters are 
chosen to bring out the fact that their effect vanishes as 
k I aH — > and B — > so long as the dynamics guarantees 
stable behavior of the parameters themselves (see giV|) . 

The Einstein equations can then be recast as a second 
order differential equation for e and constraint equations 
for the other metric variables. Since e itself contains sec- 
ond derivatives of the fundamental metric perturbations 
$ and ty, the equations are implicitly fourth order. It will 
be convenient to separate out two linear combinations of 
the underlying metric fluctuation 

$_ = -($-*), S = 2$ + *, (26) 



and a reduced mass scale or rescaling of G 



/i 2 (lna) = 



(27) 



l + f R Qna) " 

In terms of these variables, the Oi component of the Ein 
stein equations becomes a dynamical equation for e 



B{e' + Gxe) = ^G 2 S 



1 
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where 

E> E" AE" + E 
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(29) 



The metric fluctuations S and <f>_ act as sources to e 
which then feed back into their evolution weighted by 
(Bk/aH) 2 . To complete this system, the dynamics of 6 
are supplied by the derivative of its definition Eq. 
combined with Eq. (|25[) 



3 E' B' 

~ Ye + ~b 



= 2*e. 
E' 



(30) 



Eqs. (|28|) and (|30|) can also be combined to eliminate 9 
leaving a second-order differential equation for e. This 
combined relation may alternately be derived directly 
from the trace of the ij component of the Einstein equa- 
tions. 

As in general relativity, the remaining Einstein equa- 
tions become constraint equations given the dynamical 
variables A, Hq, e, 9. The 00 equation may be expressed 
as the modified Poisson equation, 



2$_ 



BE' 



2 E AE' + E 



E ' jjiS + ZB^^^A, (31) 



k 2 



where A is the density perturbation in the comoving 
gauge 



A = S - 3Hq , 
and the trace-free ij component becomes 



(32) 



$ + * 



2 E v H ' \aH 



(33) 
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E' 
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-(S + 3Be) 



--B9 

2 E 



SAE' 

Note that as B — > and £l — > fi these constraint equations 
become the usual Poisson and anisotropy equations. In 
particular for B = 0, the dark energy-closure relation 
<f> = — ^ is recovered. 

Finally, the conservation laws provide the dynamics for 
the matter fluctuations and are unmodified by f(R) 



A' 



k 



Hq - 3C , 



Hq' = = --(5-4$, 



(34) 



(35) 



The impact of the modification to gravity comes from the 
metric evolution. Eq. (|24[) implies 



S H \aH 



B9. 



(36) 



In fact directly integrating Eq. (|36[) and checking for 
consistency between Hq defined through Eq. (|19[) and 
Eq. ([55]) tests the numerical accuracy of solutions. 

Along with initial conditions for each of the fluctu- 
ations, these equations provide a complete and exact 
description of scalar linear perturbation theory in f(R) 
gravity for a matter-only universe. 
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IV. STABILITY AT HIGH CURVATURE 

The fourth-order nature of the linear perturbation 
equations derived in the previous section raises the ques- 
tion of stability in the high-curvature limit to general 
relativity [H, [53, [SH, H2] . Strongly unstable metric fluc- 
tuations can create order unity effects that invalidate the 
background expansion history. 

The key equations for stability are (|28|) and ([30| which 
describe the evolution of the deviation parameters. Con- 
sider the high rcdshift limit of high curvature where 
\B\ — *■ and wavelengths of interest are well outside the 
horizon k/aH <C 1. In this limit the evolution equations 
simplify to 



7 -+4 PB )e' + ^e=±F(<S>^,S,Hq), (37) 



where .F(<i>_, S, Hq) is the source function for the de- 
viation e and recall ps is the growth index of B from 
Eq. (TT8"|) . The details of F are not important for the sta- 
bility analysis other than that it provides a source that 
is of order the perturbation parameters that are its argu- 
ments. Under the assumption that the general relativis- 
tic solution is stably recovered in this limit, it acts as an 
external source to the deviations. The stability question 
can be phrased as whether e remains self-consistently of 
order these sources or grows and prevents the recovery of 
the solutions. 

The stability equation ([37]) has the peculiar feature 
that the frequency squared 2/B diverges as \B\ — > in- 
dependently of fc, resembling a divergent real or imagi- 
nary mass term. Evolution of e can occur on a time scale 
much shorter than the expansion time. If B < 0, e is 
highly unstable and deviations will grow exponentially. 
If B > 0, e is highly stable and is driven to the value 
required by the source function e = F/2. This short time 
scale behavior can also be seen directly in the 4th order 
form of the Einstein equation. The trace of the ij equa- 
tion or the derivative of the Oi equation have their 4th 
order terms multiplied by the small parameter B. 

Thus despite the apparent recovery of general relativity 
in the action at high curvature R, the general-relativistic 
solutions to linear perturbation theory are not recovered 
for B < 0. In terms of f{R), B cx fun in this limit and 
hence models like [l[ 



f(R) = 



M 2 + 2 



(n > 0) . 



and ml 



f(R) = —M 2 exp(-R/XM 2 ) 



(38) 



(39) 



are included in this class of unstable models. 

The instability causes any finite patch of a universe 
that starts at high curvature to break away from the 
background solution into either a low curvature solution 
or a singularity. The low curvature R <C Gp solutions 



on the other hand are stable and in fact correspond to 
the background expansion histories studied by [H, [33[ . 
However these expansion histories have gravity modified 
throughout the matter dominated epoch and in partic- 
ular a cx t 1 / 2 . They produce phenomenology at high 
redshift that would violate constraints from the CMB. 
We will omit them from further consideration below. 



V. METRIC EVOLUTION SOLUTIONS 

In this section we discuss the numerical solutions of the 
linear perturbation equations on the stable B > branch. 
To expose the underlying features of the solutions, we ex- 
amine the relevant limiting cases below. We begin with 
the initial epoch where \B\ <C 1 and the fluctuations are 
super horizon sized k/aH <C 1. We then examine large 
scale or "superhorizon" modes where B 1 / 2 k/aH <C 1 
whose evolution is completely determined by the back- 
ground expansion history and the form of f{R). Finally 
we track the evolution of small scale or "subhorizon" 
modes until B 1 / 2 k/aH ^> 1 where their evolution reaches 
the simple form implied by quasistatic equilibrium. 



A. Initial Conditions 

On the stable B > branch, we can set the initial 
conditions when B <C 1 and the mode is superhorizon 
sized, k/aH -C 1. In this case, the initial conditions 
for the normal fluctuation parameters follow the general- 
relativistic expectation 





3, 


$i = 






*i = 




A, = 


2 / k 




3 \aH 


H Qi = 





(40) 

where £;=const. is the initial comoving curvature. These 
relations also imply $_ = $ and 5 = $ with vanishing 
first derivatives initially. 

Detailed balance gives the deviation parameters as 

\pb<S>i, 



9' 



PB 0i , 



(41) 



and the high frequency term in their evolution equations 
ensures that they stay locked to these relations until B 
becomes non-negligible. 



B. Superhorizon Evolution 

Given that 9 and e are locked to the initial values of 
Eq. (|4"Tj) when B <C 1, their definitions in Eq. 
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imply that they have negligible effect on the evolution of 
the metric fluctuations <&, This remains true even as 
the mode evolves into the B ~ 1 regime if k/aH <C 1. In 
particular, the anisotropy relation of Eq. becomes 



$ + * = BH'q 



(42) 



and closes the general relation for superhorizon metric 
fluctuations Eq. (|2"3"]l: 



„ / J?" B' 

$ + 1 1 

H' 1 — B 



B^)<I>' (4:-!) 



iT _ iT B' 
~H ~ ~H T + 1 - B 



$ = 0, (fc = 0). 



The evolution of <& is completely determined by the back- 
ground evolution and the specification of f(R) ■ Formally, 
this solution also applies to the unstable branch B < 
at k = but is only valid at finite k for large \B\. The 
point at which B — 1 is a regular singular point for typi- 
cal B(\na) and so $ evolves smoothly through it. $ will 
grow on the expansion time scale if 



El 



B' 
l-B 



<0. 



(44) 



Growth is typically a transient phenomenon at the onset 
of acceleration given that the presence of matter makes 
H 1 /H — H" /H' positive. For example, if the expansion 
rate approaches the de Sitter case of a constant in the 
future H' /H — * and the solution to Eq. ([44]) becomes 



$ = 



Ci 



9i 

a 



da{l-B)H', (H'/H^O), (45) 



where C\ and C% are constants. This implies decaying 
solutions unless BH' grows. H' decays and, since in de 
Sitter R — > const, B should asymptote to a constant. 
Eq. (j4"5|) then implies <I> oc 1/a. Note that this stability 
analysis differs from treatments that take a pure de Sitter 
expansion with no matter since that assumption forces a 
closure relation of $ = * (c.f. [H, HH). 

An example of the superhorizon evolution of the metric 
is shown in Fig.[2](fc = curves) for a model with Bo = 1. 
While $ is monotonically smaller than the B = ACDM 
prediction, <f>_ is monotonically larger due to the closure 
relation between $ and of Eq. 



Subhorizon Evolution 



For subhorizon scales where k/aH ^§> 1, Eqs. (|28p and 
pop form an oscillator equation whose frequency scales 
as k/aH. Therefore the amplitude of e is driven to zero 
when compared with When combined with the Pois- 
son and anisotropy equation, this requires [3l| 



lim S 

Bk/aH^oo 



0, *^-2$, $. 



-4> 



(46) 
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FIG. 2: Evolution of metric fluctuations $ (upper panel) and <t>_ 
(lower panel) for Bo = 1 and a ACDM expansion history. The 
different closure relations on super and sub- horizon scales for 
Eqs. H42B and 1461 1, lead to qualitatively different evolution for the 
two limits with a transition region in between. <f>_, which controls 
effects in the CMB and enters directly into the Poisson equation, 
has a scale-dependent growth that makes it increasingly larger than 
the ACDM prediction at high k. Results for other values of Bo can 
be scaled from this figure by noting that the transition occurs when 
k/aH « B -1 / 2 . 



The Poisson equation then takes the simple form 



1 



il 2 a 2 pA 



(47) 



and the conservation laws become 



A' 



k 
~aH 



Hq, 



Hq' = 
q 3 ' 



(48) 



where we have dropped temporal derivatives when com- 
pared with spatial gradients where appropriate. This sys- 
tem describes a scale free evolution for $_ or A. The 
transition between these two scale-free regimes occurs 
when (k/aH) ~ B~ x l 2 . This scale- and time-dependent 
transition leads to a scale-dependent growth rate. Un- 
like for $, has monotonically enhanced power as k 
increases on the B > branch. Because of the time de- 
pendence of the transition, the total growth to z — 
continues to increase with k even for k/aH 3> B~ x l 2 . 

In Fig. [21 we show the full numerical solution from the 
initial conditions through the super- to the sub-horizon 
evolution for a few representative modes. 
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VI. POWER-SPECTRA OBSERVABLES 

The scale dependences of the linear growth rate of met- 
ric and density perturbations change predictions for cos- 
mological power spectra in the linear regime. Let us make 
the usual assumption that the initial spectrum of fluctu- 
ations in the comoving curvature is given by a power 
law. For a starting epoch during matter domination, 
this power law is modified by the usual matter-radiation 
transfer function T(k) 



k 3 P a 
2tt 2 



(49) 



where 5^ is the rms amplitude at the normalization scale 
k n . 

The modifications to the CMB depend on the scale- 
dependent change in the potential growth rate 



G(a,k) = 



$_(a,jfe) 
$_(ai,fc) 



(50) 



through the Integrated Sachs- Wolfe (ISW) effect. This 
effect comes from the differential redshift that CMB pho- 
tons suffer as they transit the evolving potential. It con- 
tributes to the angular power spectrum of temperature 
anisotropics as 



C 



47T 



k [hm 25 2tt2 



where 



l](k) = 2 J dz 



dG 

dz 



Ji(kD). 



(51) 



(52) 



Here D = J dz/H is the comoving distance out to red- 
shift z. 

In Fig. [31 we show the quadrupole power as a func- 
tion of Bq contributed by the ISW effect as well as the 
total quadrupole. Power in the quadrupole arises near 
scales of k/Ho ~ 10 and so the weak evolution of <!>_ 
shown in Fig. [2] reduces it near Bq ~ 1. In fact there 
is a minimum around Bq ss 3/2, where the ISW effect 
is a negligible contributor to the power, and a substan- 
tial reduction for 0.2 < Bq < 2.5. Further reduction of 
large-scale power can be achieved by changing the ini- 
tial power spectrum to simultaneously suppress horizon 
scale power in the Sachs- Wolfe effect from recombina- 
tion. Hence these models provide the opportunity to 
bring the predicted ensemble-averaged quadrupole closer 
to the measurements on our sky [25j . Models with Bq > 3 
produce an excess of large-angle anisotropy and exacer- 
bate the tension with the data. Note however that due to 
sample variance, changes in the likelihood will be small. 
We will address constraints on the models in a separate 
work. 

We show the full spectrum of temperature anisotropy 
Cf T in Fig.|3]for a few representative values of Bq. Given 



that the changes to the power spectrum occur mainly at 
the lowest multipoles, WMAP constraints on the ampli- 
tude of the peaks can be directly translated into a nor- 
malization of the power spectrum on scales corresponding 
to the acoustic peaks. For the ACDM expansion history 
of w = —1, f^DE = 1 — = 0.76, h = 0.73 the normal- 
ization from WMAP is 5^ = 4.52 x 10~ 5 for an optical 
depth to reionization of r = 0.092. We further take a tilt 
of n = 0.958, fl b h 2 = 0.0223. 




FIG. 3: CMB quadrupole power 6C2/2n contributed by the modi- 
fied ISW effect (dashed curve) and total (solid curve) as a function 
of Bq in the ACDM expansion history. For reference, the ACDM 
total quadrupole is also shown (horizontal line). The change in 
the growth of the potential causes a near nulling of the ISW ef- 
fect at Bo 3/2 and a substantial reduction of power between 
0.2 < B < 2.5. 




FIG. 4: CMB angular power spectra for the ACDM expansion 
history for Bo = (ACDM), 1/2, 3/2. Power in the low multipoles 
is lowered by the reduction in the ISW effect. Power at the high 
multipoles of the acoustic peaks is left unchanged. 



The WMAP normalization then allows us to predict 
the matter power spectrum today. Let us define the den- 
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sity growth rate 



Dc(a, k) = 



A(a,k) 
A(ai,k)' 



(53) 



such that Dq = a before f(R) effects become important. 
In f(R) models the potential and density growth rates 
Eq. (fST)|) and ([53")) can differ non-trivially due to the time 
dependent (1+fn) rescaling of G in the Poisson equation. 
The linear power spectrum then becomes 



2^ 



2 PL(k,a) 



25 



D 2 G (a,k) 



k 3 P, 



SI 2 H 4 



2vr 2 



(54) 
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FIG. 5: Linear matter power spectrum for several values of Bo 
in the ACDM expansion history. The change in the amplitude 
of the power at high Bo > 0.1 is nearly degenerate with galaxy 
bias. Smaller values of 0.001 < Bo < 0.1 change the shape of 
the linear power spectrum at a potentially observable level. 
All spectra are normalized to the WMAP anisotropy from 
recombination. 



We show Pl(/c) for several choices of Bq in Fig. [5] De- 
spite the large change in amplitude at high k, the high 
Bq models cannot be automatically ruled out by galaxy 
clustering data since the nearly multiplicative shift can be 
mimicked by galaxy bias. Likewise, non-linear measure- 
ments of the mass power spectrum through the cluster 
abundance, Lyman-a forest, and cosmic shear also can- 
not be straightforwardly applied. As the local curvature 
exceeds the background curvature in collapsed dark mat- 
ter halos one would expect the gravitational dynamics to 
return to Newtonian. For this reason, our predictions are 
restricted to the linear regime at k < O.lh Mpc -1 . We 
intend to explore these issues further in a future work. 

Finally the cross correlation between the ISW effect 
and the angular power spectra of galaxies is markedly 
different in these models and potentially excludes large 
Bq solutions. The angular power spectra of galaxies is 
given in the linear regime by 



fc 4 k 3 P r 



25 n 



2 ZT4 

-"o 



2tt 2 



(55) 



where 

Jf , (A) = J ' dzD G {a,k)n ] {z)b {z) ] i{kD) , (56) 

nj(z) is the galaxy redshift distribution normalized to 
J dzrij — 1, and bj(z) is the galaxy bias. 

The cross correlation between the CMB ISW effect and 
galaxies becomes 



Cf = ^j d ^(k)l]{k) 



6 



k 3 Pr 



(57) 



The correlation coefficient between the total temperature 
anisotropy and the galaxies is given by 



Ri 



(jTT(j9i9i 



(58) 



and is independent of the galaxy bias if it is slowly vary- 
ing with redshift. We have neglected galaxy magnifica- 
tion bias which leads to an additional source of correla- 
tion. 

For definiteness, we assume that the galaxy sets come 
from a net galaxy distribution of 



n g (z) oc z 2 e -t*A-5) 2 , 



(59) 



which is further partitioned by photometric redshift into 
several galaxy samples, 



n {z) oc n g (z) 



crfc 



Z 3-l 



— Z 



V2a(z 



erfc 



( Z 3~ Z 

\V2<t(z) 



(60) 



where erfc is the complementary error function and 
cr(z) = 0.03(l-|-z) reflects the effect of photometric red- 
shift scatter. 

Fig. [6] shows the correlation coefficient for several val- 
ues of Bq. We take two redshift bins from Eq. (foT))) par- 
titioned by Zj — 0, 0.4, 0.8 to achieve effective redshifts 
of z = 0.2 and 0.6. Current observations constrain the 
correlation near / ~ 20 corresponding to scales which 
are an order of magnitude smaller than those contribut- 
ing to the ISW quadrupole. Between 1/2 < B < 1 
the galaxy-ISW correlation is substantially reduced. For 
Bq > 3/2, galaxies are in fact anti-correlated with the 
CMB since $_ grows on the relevant scales. A loose 
bound from the observed correlation would therefore be 
Bo ^ 3/2 at th e sig nificance levels of the reported detec- 
tions (e.g. [H, HJIH H^, H3, EH) but we expect more 
detailed modeling to yield better constraints in the fu- 
ture. It is likely that a significant reduction of the large 
angle anisotropy from this mechanism could be excluded 
unless other sources, such as magnification bias, can gen- 
erate the observed positive correlation. 



VII. DISCUSSION 

We have studied the evolution of linear cosmological 
perturbations in f(R) models for modified gravity in the 
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FIG. 6: Cross correlation coefficient between the CMB and 
galaxies in the ACDM expansion history. Shown are two rep- 
resentative redshift bins centered around z = 0.2 and 0.6 with 
Bo = (ACDM), 1/2, 1, 3/2. The cross correlation is sub- 
stantially reduced for 1/2 < Bo < 1 and becomes negative for 
Bo > 3/2. 



with CMB measurements. 

For the stable B > branch, f(R) models predict a 
rich set of linear phenomena that can be used to test 
such deviations from general relativity. For example, 
large B ~ 1 models lower the large-angle anisotropy 
of the CMB and may be useful for explaining the low 
quadrupole observed on our sky. They also predict qual- 
itatively different correlations between the CMB and 
galaxy surveys which may provide the best upper limit 
on the deviations currently available. Smaller deviations 
in B are observable at smaller scales through changes to 
the shape of the linear power spectrum. In the limit that 
B — ► and the expansion history is given by ACDM, lin- 
ear perturbations in f(R) models approach the general 
relativistic predictions exactly. We intend to examine 
constraints on f(R) models in a future work. 

More generally, this class of phenomenological f(R) 
models provides insight on the types of deviations that 
might be expected from alternate metric theories of grav- 
ity in the linear regime. Conservation of the matter 
stress-energy tensor severely restricts the form of allowed 
deviations on both super- and sub- horizon scales fl6ll27j|. 
Even if these f(R) models prove not to be viable as a 
complete alternate theory of gravity that includes solar- 
system tests, they may serve as the basis for a "parame- 
terized post-Friedmann" description of linear phenomena 
that parallels the parameterized post-Newtonian descrip- 
tion of small-scale tests. 



physical — or Jordan — frame. Here the gravitational dy- 
namics are fourth order and the matter is minimally cou- 
pled and separately conserved. For models that recover 
the Einstein-Hilbert action at high curvature R, we find 
that for each expansion history specified by H (In a) there 
exists two branches of f(R) solutions that are parame- 
terized by B oc /rr, the second derivative of f(R). For 
B < 0, which includes most models previously considered 
P, [3l| , there is a short-timescale instability that prevents 
recovery of the general-relativistic expectations at high 
curvature that is important for maintaining agreement 
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